Three different methods to determine the concentration 
profile for a tubular reactor with axial dispersion
Author's Data: Housam Binous
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Acknowledgement :

This problem was presented in the excellent book by W. F. Ramirez, Computational Methods for Process Simulation, Second Edition, 1997, Butterworth Heinemann, Oxford.

The shooting method has been developped using Mathematica by Professor Brian Higgins (Chemical Engineering Department, University of California, Davis) is available at : http://www.higgins.ucdavis.edu/chemmath.php .

Off[General::"spell1"]

Guessing the value of CA[y=0] using three methods

First Method inspired by Professor Fred Ramirez

CA0 = 0.5;

delta = 10;

i = 1;

While[
      delta > 0.0001, {sol = 
          NDSolve[{D[CA[y], {y, 1}] ==  -x1[y], 
              D[x1[y], {y, 1}] ==  -100 x1[y] - 200 CA[y], 
              x1[0] == 0, CA[0] == CA0, 
              D[f1[y], {y, 1}] ==  -f2[y], 
              D[f2[y], {y, 1}] ==  -200 f1[y] - 100 f2[y], 
              f1[0] == 1, f2[0] == 0}, {CA[y], x1[y], f1[y], f2[y]}, {y, 0, 
              1}], 
delta = ((1 - CA[y] + 
                        0.01 x1[y])/(-f1[y] + 0.01 f2[y]) /. 
                  sol /. y -> 1)[[1]], CA0 += -delta, 
        i = i + 1, Print[delta, " ", i]}]

CA0

Second Method inspired by Professor Brian Higgins (U C Davis)

CA0 =.

SOLD[CA0_] := 
      NDSolve[{D[CA[y], {y, 1}] ==  -x1[y], 
          D[x1[y], {y, 1}] ==  -100 x1[y] - 200 CA[y], 
          x1[0] == 0, CA[0] == CA0}, {CA[y], x1[y]}, {y, 0, 
          1}]

yend[CA0_?NumericQ] := 
      1 - (CA[y] /. SOLD[CA0] /. 
            y -> 1) + (0.01 x1[y] /. SOLD[CA0] /. 
            y -> 1)

FindRoot[yend[CA0] ==  0, {CA0, 0.5}, 
      MaxIterations ->  1000]

Plotting the concentration profile

First[CA[y] /. SOLD[0.140592] /. y -> 0]

f = CA[y] /. SOLD[0.140592]

plt1 = 
      Plot[(f /. y -> 1 - x), {x, 0, 1}, 
        PlotStyle -> RGBColor[1, 0, 0]]

Third Method using NDSolve

solpde = 
      NDSolve[{D[CA[x], {x, 2}] - 100 D[CA[x], {x, 1}] - 200 CA[x] == 
            0, 
CA'[1] ==  0, 
          1 ==  CA[0] - 0.01 CA'[0]}, {CA[x]}, {x, 0, 
          1}]

plt2 = 
      Plot[Evaluate[CA[x] /. solpde], {x, 0, 1}, 
        PlotStyle -> RGBColor[0, 1, 0]]

CA[x] /. solpde /. x -> 1

Comparing concentration profiles obtained using the three methods

Show[plt1, plt2]